Search for: All records

Abstract Megathrust earthquakes are the largest on Earth, capable of causing strong ground shaking and generating tsunamis. Physical models used to understand megathrust earthquake hazard are limited by existing uncertainties about material properties and governing processes in subduction zones. A key quantity in megathrust hazard assessment is the distance between the updip and downdip rupture limits. The thermal structure of a subduction zone exerts a first‐order control on the extent of rupture. We simulate temperature for profiles of the Cascadia, Nankai and Hikurangi subduction zones using a 2D coupled kinematic‐dynamic thermal model. We then build reduced‐order models (ROMs) for temperature using the interpolated Proper Orthogonal Decomposition (iPOD). The resulting ROMs are data‐driven, model agnostic, and computationally cheap to evaluate. Using the ROMs, we can efficiently investigate the sensitivity of temperature to input parameters, physical processes, and modeling choices. We find that temperature, and by extension the potential rupture extent, is most sensitive to variability in parameters that describe shear heating on the slab interface, followed by parameters controlling the thermal structure of the incoming lithosphere and coupling between the slab and the mantle. We quantify the effect of using steady‐state versus time‐dependent models, and of uncertainty in the choice of isotherm representing the downdip rupture limit. We show that variability in input parameters translates to significant differences in estimated moment magnitude. Our analysis highlights the strong effect of variability in the apparent coefficient of friction, with previously published ranges resulting in pronounced variability in estimated rupture limit depths. 

Abstract Reconstructing fault surfaces from volumetric data is a longstanding challenge in geosciences. We present a novel 3D method based on the medial axis to transform a volumetric strain‐rate invariant field from long‐term geodynamic simulations into fault surfaces. In these geodynamic models, faults correspond to regions of locally high values of the second invariant of the strain‐rate commonly referred to as shear zones. The proposed workflow begins by normalizing the strain‐rate to define fault indicator field . An iso‐surface of a chosen value is then extracted to form an envelope around the shear zones. Using the shrinking ball algorithm (Ma et al., 2012,https://doi.org/10.1007/s00371‐011‐0594‐7), we compute the medial axis of this 3D envelope to generate a point cloud representing the geometric skeleton of the shear zones. We reconstruct fault surfaces by applying Delaunay triangulation followed by Laplacian smoothing. For models involving multiple intersecting faults, we perform a local principal component analysis (PCA) of the coordinates defining the medial axis and use the resulting eigenvectors to detect first‐order orientation variations, enabling the separation and individualization of faults. We demonstrate the generality and robustness of the method by applying it several diverse 3D geodynamic scenarios: A single strike‐slip fault, a branching strike‐slip fault in a restraining bend, a dense strike‐slip fault network, a rift system, and a subduction zone with a megathrust and a conjugate thrust fault. 

Abstract Structural inversion of rifted basins is generally associated with surface uplift and denudation of the sedimentary infill, reflecting the active contractional deformation in the crust. However, worldwide examples of inverted rifts show contrasting basin-scale subsidence and widespread sedimentation patterns during basin inversion. By conducting a series of three-dimensional coupled geodynamic and surface processes models, we investigated the dynamic controls on these subsidence anomalies during the successive stages of rifting and basin inversion, and we propose a new evolutionary model for this process. Our models show that the inherited thermo-rheological properties of the lithosphere influence the initial strain localization and subsequent migration of crustal deformation during inversion. The sense of the vertical movements (i.e., uplift or subsidence), however, is not directly linked to the underlying crustal stress patterns; rather, it reflects the balance among contraction-induced tectonic uplift, postrift thermal subsidence of the inherited lithosphere, and sediment redistribution. Based on the interplay among the competing differential vertical movements with different amplitudes and wavelengths, inversion of rifted basins may lead to the growth of intraplate orogens, or the contraction-driven localized uplift may be hindered by the thermal sag effects of the inherited shallow lithosphere-asthenosphere boundary, resulting in basin-scale subsidence. In such basins, dating the first erosional surfaces and other unconformities may not provide accurate timing for the onset of inversion. 

SUMMARY To reach Earth’s surface, magma must ascend from the hot, ductile asthenosphere through cold and brittle rock in the lithosphere. It does so via fluid-filled fractures called dykes. While the continuum mechanics of ductile asthenosphere is well established, there has been little theoretical work on the cold and brittle regime where dyking and faulting occurs. Geodynamic models use plasticity to model fault-like behaviour; plasticity also shows promise for modelling dykes. Here we build on an existing model to develop a poro-viscoelastic–viscoplastic theory for two-phase flow across the lithosphere. Our theory addresses the deficiencies of previous work by incorporating (i) a hyperbolic yield surface, (ii) a plastic potential with control of dilatancy and (iii) a viscous regularization of plastic failure. We use analytical and numerical solutions to investigate the behaviour of this theory. Through idealized models and a comparison to linear elastic fracture mechanics, we demonstrate that this behaviour includes a continuum representation of dyking. Finally, we consider a model scenario reminiscent of continental rifting and demonstrate the consequences of dyke injection into the cold, upper lithosphere: a sharp reduction in the force required to rift. 

SUMMARY Physics-based simulations provide a path to overcome the lack of observational data hampering a holistic understanding of earthquake faulting and crustal deformation across the vastly varying space–time scales governing the seismic cycle. However, simulations of sequences of earthquakes and aseismic slip (SEAS) including the complex geometries and heterogeneities of the subsurface are challenging. We present a symmetric interior penalty discontinuous Galerkin (SIPG) method to perform SEAS simulations accounting for the aforementioned challenges. Due to the discontinuous nature of the approximation, the spatial discretization natively provides a means to impose boundary and interface conditions. The method accommodates 2-D and 3-D domains, is of arbitrary order, handles subelement variations in material properties and supports isoparametric elements, that is, high-order representations of the exterior boundaries, interior material interfaces and embedded faults. We provide an open-source reference implementation, Tandem, that utilizes highly efficient kernels for evaluating the SIPG linear and bilinear forms, is inherently parallel and well suited to perform high-resolution simulations on large-scale distributed memory architectures. Additional flexibility and efficiency is provided by optionally defining the displacement evaluation via a discrete Green’s function approach, exploiting advantages of both the boundary integral and volumetric methods. The optional discrete Green’s functions are evaluated once in a pre-computation stage using algorithmically optimal and scalable sparse parallel solvers and pre-conditioners. We illustrate the characteristics of the SIPG formulation via an extensive suite of verification problems (analytic, manufactured and code comparison) for elastostatic and quasi-dynamic problems. Our verification suite demonstrates that high-order convergence of the discrete solution can be achieved in space and time and highlights the benefits of using a high-order representation of the displacement, material properties and geometries. We apply Tandem to realistic demonstration models consisting of a 2-D SEAS multifault scenario on a shallowly dipping normal fault with four curved splay faults, and a 3-D intersecting multifault scenario of elastostatic instantaneous displacement of the 2019 Ridgecrest, CA, earthquake sequence. We exploit the curvilinear geometry representation in both application examples and elucidate the importance of accurate stress (or displacement gradient) representation on-fault. This study entails several methodological novelties. We derive a sharp bound on the smallest value of the SIPG penalty ensuring stability for isotropic, elastic materials; define a new flux to incorporate embedded faults in a standard SIPG scheme; employ a hybrid multilevel pre-conditioner for the discrete elasticity problem; and demonstrate that curvilinear elements are specifically beneficial for volumetric SEAS simulations. We show that our method can be applied for solving interesting geophysical problems using massively parallel computing. Finally, this is the first time a discontinuous Galerkin method is published for the numerical simulations of SEAS, opening new avenues to pursue extreme scale 3-D SEAS simulations in the future. 

Abstract. Modelling the pressure in the Earth's interior is a common problem in Earth sciences. In this study we propose a method based on the conservation of the momentum of a fluid by using a hydrostatic scenario or a uniformly moving fluid to approximate the pressure. This results in a partial differential equation (PDE) that can be solved using classical numerical methods. In hydrostatic cases, the computed pressure is the lithostatic pressure. In non-hydrostatic cases, we show that this PDE-based approach better approximates the total pressure than the classical 1D depth-integrated approach. To illustrate the performance of this PDE-based formulation we present several hydrostatic and non-hydrostatic 2D models in which we compute the lithostatic pressure or an approximation of the total pressure, respectively. Moreover, we also present a 3D rift model that uses that approximated pressure as a time-dependent boundary condition to simulate far-field normal stresses. This model shows a high degree of non-cylindrical deformation, resulting from the stress boundary condition, that is accommodated by strike-slip shear zones. We compare the result of this numerical model with a traditional rift model employing free-slip boundary conditions to demonstrate the first-order implications of considering “open” boundary conditions in 3D thermo-mechanical rift models. 

Abstract In traditional modeling approaches, earthquakes are often depicted as displacement discontinuities across zero‐thickness surfaces embedded within a linear elastodynamic continuum. This simplification, however, overlooks the intricate nature of natural fault zones and may fail to capture key physical phenomena integral to fault processes. Here, we propose a diffuse interface description for dynamic earthquake rupture modeling to address these limitations and gain deeper insight into fault zones' multifaceted volumetric failure patterns, mechanics, and seismicity. Our model leverages a steady‐state phase‐field, implying time‐independent fault zone geometry, which is defined by the contours of a signed distance function relative to a virtual fault plane. Our approach extends the classical stress glut method, adept at approximating fault‐jump conditions through inelastic alterations to stress components. We remove the sharp discontinuities typically introduced by the stress glut approach via our spatially smooth, mesh‐independent fault representation while maintaining the method's inherent logical simplicity within the well‐established spectral element method framework. We verify our approach using 2D numerical experiments in an open‐source spectral element implementation, examining both a kinematically driven Kostrov‐like crack and spontaneous dynamic rupture in diffuse fault zones. The capabilities of our methodology are showcased through mesh‐independent planar and curved fault zone geometries. Moreover, we highlight that our phase‐field‐based diffuse rupture dynamics models contain fundamental variations within the fault zone. Dynamic stresses intertwined with a volumetrically applied friction law give rise to oblique plastic shear and fault reactivation, markedly impacting rupture front dynamics and seismic wave radiation. Our results encourage future applications of phase‐field‐based earthquake modeling. 

Abstract Physics‐based simulations of earthquake ground motion are useful to complement recorded ground motions. However, the computational expense of performing numerical simulations hinders their applicability to tasks that require real‐time solutions or ensembles of solutions for different earthquake sources. To enable rapid physics‐based solutions, we present a reduced‐order modeling approach based on interpolated proper orthogonal decomposition (POD) to predict peak ground velocities (PGVs). As a demonstrator, we consider PGVs from regional 3D wave propagation simulations at the location of the 2008M_W5.4 Chino Hills earthquake using double‐couple sources with varying depth and focal mechanisms. These simulations resolve frequencies ≤1.0 Hz and include topography, viscoelastic attenuation, and S‐wave speeds ≥500 m/s. We evaluate the accuracy of the interpolated POD reduced‐order model (ROM) as a function of the approximation method. Comparing the radial basis function (RBF), multilayer perceptron neural network, random forest, andk‐nearest neighbor, we find that the RBF interpolation gives the lowest error (≈0.1 cm/s) when tested against an independent data set. We also find that evaluating the ROM is 10⁷–10⁸times faster than the wave propagation simulations. We use the ROM to generate PGV maps for 1 million different focal mechanisms, in which we identify potentially damaging ground motions and quantify correlations between focal mechanism, depth, and accuracy of the predicted PGV. Our results demonstrate that the ROM can rapidly and accurately approximate the PGV from wave propagation simulations with variable source properties, topography, and complex subsurface structure.